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1.0  PROJECT  STATUS 


The  work  on  this  project  started  on  May  20,  1997.  The  project  is  currently  on  schedule. 

2.0  DESCRIPTION  OF  PROGRESS 

The  work  performed  during  this  quarter  is  listed  below: 

(i)  An  Implicit  Multi-dimensional  Model  for  Electroosmosis  has  been  developed  and 
successfully  integrated  with  CFD-ACE+  code. 

(ii)  Extension  of  the  above-mentioned  model  to  simulate  electrophoretic  flow  is  currently 
underway. 

The  following  sections  will  provide  detailed  discussions  on  the  activities  performed  during  this 
reporting  period. 

3.0  AN  IMPLICIT  MULTIDIMENSIONAL  FINITE-VOLUME  MODEL  FOR 
ELECTROOSMOTIC  FLOW 

3.1  Introduction 

As  mentioned  in  the  previous  quarterly  reports,  a  one-dimensional  implicit  finite-volume  method 
to  model  electrophoresis  phenomena  has  been  developed.  In  the  present  reporting  period,  this 
model  has  been  extended  for  multi-dimensional  geometries  to  simulate  electroosmotic  flow  in 
microdevices.  Electrophoresis  model  in  CFD-ACE+  will  be  incorporated  during  the  next 
quarterly  period. 

3.2  Governing  Equations 

Mass  transport  of  a  charged  specie  due  to  electroosmosis  is  governed  by  three  major 
mechanisms,  such  as  diffusion,  convection  and  electric  field.  To  model  these  mechanisms,  the 
equations  governing  the  convection  and  diffusion  of  charged  species  in  an  electric  field  must  be 
solved.  The  governing  equations  are: 

Continuity  Equation 


v-v  =  o 


(1) 


Momentum  Equation 
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where  V  is  the  velocity  vector,  p  is  the  density  of  the  buffer  solution,  p  is  pressure,  |i  is  the 
dynamic  viscosity  of  the  buffer,  pe  is  the  electric  charge  density,  and  E  is  the  electric  field 
intensity.  The  electric  field  is  related  to  the  electric  potential  <E>  via 

E  =  -V<j)  (3) 

The  electric  potential  is  governed  by 


V2$  =  -(pe/£)  (4) 

where  s  is  the  electrical  permittivity  of  the  buffer  solution.  Substituting  Eqs.  (3)  and  (4)  into  (2), 
the  momentum  equation  can  be  written  as: 


fr+(V'v)v 


=  -  Vp  +  pv2v  +  e  jv  2o  •  vo] 


(5) 


In  the  above  formulation  thermal  convection  due  to  Joule  heating,  and  chemical  reactions  within 
buffer  solution  have  been  neglected.  In  the  electrolyte  solution,  positively  and  negatively 
charged  ions  resulting  from  the  reactions  of  dissociation  interact  with  the  electric  field  that  can 
be  calculated  from  Eq.  (4).  On  a  macroscopic  scale,  the  solution  is  electrically  neutral.  The 
electro-neutrality  condition  is  violated  only  on  scales  comparable  to  Debye  radius  around  each 
ion.  This  radius  is  usually  of  the  order  of  several  nanometers.  Therefore,  for  macroscopic 
consideration  electro-neutrality  condition  can  be  used.  However,  special  attention  should  be 
given  to  regions  very  close  to  the  channel  wall,  where  a  electrical  double  layer  would  exist  and 
same  charge  ions  concentrate  in  close  proximity  resulting  in  a  potential  called  %”  potential.  An 
approach  used  by  Patankar  and  Hu  ( 1 998)  is  used  to  solve  Eq.  (4). 


In  electroosmotic  flow,  when  an  external  electric  field  is  applied,  the  charged  species  lying  next 
to  the  wall  will  start  moving.  Thus  the  magnitude  of  the  ^  potential  and  applied  electric  field 
determine  the  speed  of  motion.  Typically  when  the  Debye  thickness  is  small  and  charge  on  the 
wall  is  not  very  significant,  the  distribution  of  £  potential  is  affected  very  little  by  the  applied 
electric  field.  Also,  when  the  inertial  terms  are  small,  i.e.,  at  very  low  flow  velocity,  the  fluid 
motion  has  very  little  influence  on  charge  distribution.  Under  these  conditions,  the  total  potential 
can  be  decomposed  into  two  parts,  <j),  the  potential  due  to  the  applied  electric  field  and  \| /, 
potential  due  to  charge  on  the  walls: 


®  =  <t>  +  V  (6) 

Consequently,  Eq.  (4)  can  be  written  as  two  separate  equations 

V2(J)  =  0  (7) 

and 

vV  =  -(pe/e)  (8) 
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The  electric  field  that  generates  the  force  on  the  charged  species  causing  electroosmotic  flow  is 
given  by 


E  =  -  V  (<j>  +  \(/)  (9) 

Following  Patankar  and  Hu  (1998)  and  adapted  from  original  works  of  Grossman  and  Colburn 
(1992),  the  charge  density  can  be  expressed  as 

pe=-£K2\j/  (10) 


where  K  is  inverse  of  Debye  thickness.  The  governing  equations  can  be  rewritten  as: 


dV 

5t 


+  (V-V)V 


=  -  Vp  +  |J.V2V  +  £K2\j/V  <J) 


(11) 


V2(j)  =  0 


(12) 


V2\|/  =  k2\)/ 


(13) 


3.3  Numerical  Scheme 

The  governing  equations  (Eqs.  (1),  (11),  (12)  and  (13))  are  first  discretized  in  space  and  time 
using  a  finite-volume  based  methodology.  An  important  motivation  for  choosing  this 
methodology  is  that  it  has  excellent  conservation  property  that  is  satisfied  even  at  the  numerical 
level.  This  ensures  that  numerical  scheme  does  not  introduce  any  spurious  source  and  sink 
terms.  The  spatial  region  is  subdivided  into  individual  non-overlapping  computational  cells 
whose  union  will  reproduce  the  entire  computational  domain.  This  discretization  is  achieved 
using  CFD-GEOM,  the  geometric-modeling  and  grid  generation  software  developed  and 
marketed  by  CFD  Research  Corporation.  Data  transfer,  storage  and  retrieval  are  accomplished 
using  CFD-DTF  utility.  The  present  model  is  also  capable  of  handling  arbitrary  polyhedral  grids 
such  as  those  commonly  occur  in  unstructured  grids. 

The  finite  volume  discretization  transforms  the  partial  differential  equations  into  a  set  of 
algebraic  equations.  The  discretized  equations  are  assembled  into  matrices  and  solved  using 
iterative  linear  solvers.  The  solution  methodology  used  in  CFD-ACE+  is  based  on  SIMPLEC 
(Semi-Implicit  Pressure  Linked  Equations,  Consistent)  algorithm,  first  introduced  by  Van 
Doormal  and  Raithby  (1984).  This  is  derived  from  SIMPLE  algorithm,  described  in  detail  by 
Patankar  (1980). 

3.4  Geometry 

The  geometry  of  the  device  used  in  the  present  study  is  similar  to  one  used  by  Fan  and  Harrison 
(1994)  and  explained  in  detail  by  Patankar  and  Hu  (1998).  This  is  shown  in  Figure  1.  The 
device  is  shaped  like  a  cross-channel  with  the  intersection  of  four  reservoirs  numbered  1  through 
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4.  The  injection  of  the  sample  solution  is  affected  by  applying  a  voltage  between  reservoirs  1 
and  2.  This  will  drive  the  sample  from  reservoir  1  to  2  (for  <j)i>  (J>2  and  C,  <  0),  across  the 
intersection,  thus  creating  a  sample  plug  in  the  path  between  reservoirs  3  and  4.  By  applying 
suitable  electric  potential  between  reservoirs  3  and  4,  the  plug  will  be  separated  from  the  cross¬ 
stream  and  start  flowing  in  the  side  channel  connecting  reservoirs  3  and  4.  Sample  separation 
does  not  occur  during  injection  process  between  reservoirs  1  and  2  since  the  sample  flows 
continuously.  The  present  study  is  focussed  on  the  flow  patterns  developed  at  the  intersection  of 
cross  channels  under  different  applied  electric  field  and  Reynolds  number  of  the  flow.  The 
scalar  transport  equation  solver  for  sample  concentration  will  be  implemented  in  the  next  quarter. 
Since  the  flow  variations  and  magnitude  of  velocity  are  negligibly  small  in  the  z  direction,  the 
cross-channel  is  modeled  as  a  two-dimensional  geometry. 


3 


4 


Figure  1 .  Cross-Channel  Device 


Boundary  Conditions 


The  boundary  conditions  for  the  potential  equations  are  the  external  applied  voltage  and  ^ 
potential  at  the  wall.  A  constant  potential  ((pi)  is  applied  at  the  reservoir  1  and  reservoir  2  is 
grounded.  The  potential  applied  at  reservoirs  3  and  4  is  changed  to  effect  electroosmotic  flow 
and  sample  plug  injection.  The  walls  of  the  channel  are  assumed  to  be  electrically  insulated.  A 
constant  C,  potential  (space  charge  potential)  is  applied  at  the  wall.  The  driving  force  for  the  fluid 
motion  is  generated  by  the  electric  field. 


3.6  Results  and  Discussions 

The  buffer  solution  is  assumed  to  be  very  dilute: 

e  =  6.95  x  10‘10  C2/J.m 


£  =  100  Mv 

=  10‘3  N  s/m2 
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p  =  1000  Kg/m3 

4>,  =  150  V 

w  =30  jim 


where  w  is  the  width  of  the  channel.  Reynolds  number  (Re)  is  defined  as: 


Re= 


f-1 

UwL, 

(13) 


where  L  is  the  distance  between  reservoirs  1  and  2.  The  main  parameters  of  the  problem  are 
Reynolds  number  Re,  inverse  of  Debye  thickness  K,  potential  applied  at  reservoirs  3  (§3)  and  4 
(<t>4),  and  the  ratio  L/w.  The  effects  of  Reynolds  number  and  potential  at  side  reservoirs  (<|)3  and 
<f>4)  will  be  analyzed.  A  constant  potential  of  -100  mV  is  applied  at  the  walls.  This  implies  that 
the  wall  is  negatively  charged  and  will  form  a  thin  layer  of  positively  charged  ions  next  to  wall. 
Application  of  positive  electric  field  will  induce  a  fluid  velocity  in  the  positive  direction. 

3.6.1  External  Potential  Field 

Figure  2  shows  the  external  potential  field  under  the  application  of  different  potentials  at  the  side 
reservoirs  3  and  4.  Reservoir  1  is  always  kept  at  150  V  and  reservoir  2  is  always  grounded. 
When  the  potential  at  the  side  reservoirs  is  75  V,  it  is  observed  that  by  symmetry,  potential  at  the 
center  of  the  channel  is  also  75  V.  As  a  result,  there  is  very  small  variation  of  potential  in  the 
side  channel  and  the  electric  field  is  almost  zero.  There  is  no  net  flow  in  these  channels. 
Velocity  vector  for  this  case  is  shown  in  Figure  3  (a).  As  expected,  the  velocity  profile  is  flat 
(plug  flow). 


(a)  (b)  (c) 


Figure  2.  Contours  of  External  Potential  Field  Inside  Cross-Channel  for  Re  =  1.  Potential 
at  the  Side  Reservoirs:  (a)  75  V,  (b)  37.5  V,  and  (c)  135  V 
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(a)  (b)  (c) 

Figure  3.  Velocity  Vector  Plot  Inside  Cross-Channel  for  Re  =  1.  Potential  at  the  Side 
Reservoirs:  (a)  75  V,  (b)  37.5  V,  and  (c)  135  V.  The  Phenomena  of  Electrokinetic 
Focusing  Can  be  Observed  in  Figure  (c) 

When  the  potential  at  the  side  reservoir  is  37.5  V,  the  potential  at  the  center  of  the  channel  is  less 
than  75  V.  There  is  also  a  potential  gradient  along  the  side  channel  and  there  is  net  flow  away 
from  the  center  of  the  channel.  On  the  contrary,  when  the  potential  at  the  side  reservoirs  are  kept 
at  135  V,  there  is  a  net  flow  towards  the  center  of  the  channel.  This  phenomena  is  clearly 
illustrated  in  Figures  2(b)  and  (c).  In  Figures  3(b)  and  (c),  the  corresponding  velocity  vector  plot 
is  shown. 

3.6.2  Potential  Due  to  Charge  on  the  Walls 

Figure  4(a)  illustrates  a  typical  grid  used  in  the  simulation.  As  can  be  seen,  the  grid  lines  are 
densely  packed  near  the  wall  so  that  variations  in  C,  potential  due  to  wall  charge  will  be  resolved. 
In  Figure  4(b),  the  contours  of  £  potential  is  shown.  The  potential  decays  exponentially  to  zero 
at  the  center  of  the  channel.  In  all  the  simulations  shown  in  this  report,  a  constant  value  of  1.5 
pm  is  used  as  the  Debye  thickness. 


(a)  (b) 

Figure  4.  Nonuniform  Grid  Used  to  Generate  Cross-Channel  Device  is  Shown  in  (a)  and 
Contour  Plots  of  C,  Potential  is  Shown  in  (b)  for  <(>i  =  150  V  and  <f>3  =  <|>4  =  75  V 
and  Re  =  1 
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3.6.3  Effect  of  Reynolds  Number 

To  understand  the  effect  of  Reynolds  number  on  the  flow  field,  the  side  reservoirs  are  set  to  a 
potential  that  is  half  of  the  potential  applied  at  reservoir  1  (4»i).  This  will  insure  that  there  will 
not  be  any  net  flow  in  side  channels.  Also,  the  pressure  at  the  reservoirs  is  set  to  zero.  Reynolds 
number  is  varied  by  changing  <|>i.  The  velocity  vector  plot  inside  the  injection  channel  is  shown 
in  Figures  5(a)  and  (b)  for  Re  =  0.1  and  10,  respectively.  Overall  it  is  observed  that  the  velocity 
profile  remains  flat  through  out  the  channel  except  near  the  wall,  where  sharp  gradients  are 
introduced  due  to  the  imposition  of  noslip  condition.  As  Re  increases,  velocity  of  flow  in  the 
injection  channel  increases.  Consequently,  the  leakage  of  sample  from  the  injection  channel  into 
the  separation  channel  is  reduced.  This  will  result  in  more  uniform  sample  plug. 


(a)  (b) 

Figure  5.  Effect  of  Reynolds  Number  on  Flow  Field.  The  Velocity  Vector  Plot  in  Cross- 
Channel  Device  is  Shown  for  Re  =  0.1  and  10  in  (a)  and  (b)  Respectively.  The 
Potential  at  Side  Reservoirs  is  Kept  at  Half  of  <J>j.  The  Maximum  Velocity 
Observed  are  0.044  m/s  and  0.24  m/s  for  Re  =  0.1  and  10,  Respectively 

3.6.4  Electrokinetic  Focusing 

This  is  one  way  of  obtaining  tightly  focused  stream  of  samples  at  the  channel  intersection  so  that 
a  uniformly  shape  sample  plug  may  be  injected  into  side  channels  for  separation.  This  can  be 
achieved  by  varying  potential  at  the  side  reservoirs.  For  this  study,  Reynolds  Number  is  set  to 
1.0  (<|>i  =  150  V)  and  pressure  at  all  the  reservoirs  is  set  to  zero. 


As  shown  in  Figures  2(a)  and  3(a),  when  the  potential  at  the  side  reservoir  is  set  to  75  V,  there  is 
no  net  flow  in  side  channels.  When  the  potential  is  reduced  below  75  V  (to  37.5  V),  it  is 
observed  from  Figures  2(b)  and  3(b)  that  there  exists  a  potential  gradient  and  fluid  flow  occurs 
from  the  injection  channel  to  the  separation  channel.  This  is  not  desirable  since  it  may  cause  the 
sample  to  leak  from  the  injection  channel  and  may  contaminate  the  buffer  solution.  When  the 
potential  is  set  to  135  V,  it  is  evident  from  Figures  2(c)  and  3(c)  that  the  fluid  flows  from  the 
separation  channel  into  the  injection  channels.  The  streamlines  in  the  injection  channel  are 
squeezed  at  the  intersection.  This  is  desirable  since  the  sample  is  restricted  to  remain  at  the 
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intersection.  However  too  much  squeezing  can  distort  the  shape  of  the  plug.  This  technique  of 
controlling  the  shape  of  plug  by  varying  potential  at  the  side  channels  has  been  used  in 
experiments  (Patankar  and  Hu,  1998)  and  computational  software  such  as  CFD-ACE+  can  be 
used  as  a  valuable  design  tool  to  optimize  applied  potential.  More  detailed  study  of  this 
focussing  will  be  reported  in  the  next  progress  report. 

4.0  PLANNED  ACTIVITIES 

The  following  activities  are  planned  for  the  next  quarter: 

(i)  Effort  will  be  undertaken  to  integrate  the  current  electrophoresis  model  with  CFD  ACE+ 
environment. 

(ii)  Multi-dimensional  simulation  capabilities  will  be  developed  to  study  of  electrophoresis. 

(iii)  Collaborative  work  between  CFDRC  and  Aclara  Biosciences,  InnovaTek,  Biomicro 
Systems  and  Affymetrix  will  be  pursued. 

5.0  SUMMARY  OF  MEETINGS  ATTENDED 

No  meeting  has  been  attended  during  the  reporting  period 

6.0  ACCOMPLISHMENTS 

An  implicit  finite-volume  based  CFD  model  to  simulate  electroosmotic  flow  has  been  integrated 
with  CFD-ACE+. 
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PHASE  II SBIR 


ON-CHIP  TRANSPORT  OF  BIOLOGICAL  FLUIDS 
IN  MEMS  DEVICES 

DARPA  CONTRACT  NUMBER:  MDA972-97-C-0015 


Total  Funds 

$980,000 

Total  Funds  Used 

$850,461 

Estimated  Cost  to  Complete 

$129,539 

Estimated  Percentage  of  Physical  Completion 

$87% 

